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£\1 , The standard hybrid Monte Carlo algorithm is known to simulate even flavors 

QCD only. Simulations of odd flavors QCD, however, can be also performed in the 

^vj framework of the hybrid Monte Carlo algorithm where the inverse of the fermion 

^. ' matrix is approximated by a polynomial. In this exploratory study we perform 

_j_ ' three flavors QCD simulations. We make a comparison of the hybrid Monte Carlo 

— ^1 t algorithm and the R-algorithm which also simulates odd flavors systems but has 



Jl 
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step-size errors. We find that results from our hybrid Monte Carlo algorithm are 
in agreement with those from the R-algorithm obtained at very small step-size. 
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f~} 1 Introduction 

Recent lattice QCD simulations include effects of dynamical fcrmions. The 
standard established algorithm .-to simulate dynamical QCD is the Hybrid 
Q_i' Monte Carlo (HMC) algorithm tl but it is limited to simulations of an even 

number of degenerate flavors. It would be desirable to simulate lattice QCD 
with three flavors since there exist three light quarks, i.e. u,d,s quarks, in 
the real world. Simulations with an odd number of flavors have been per- 
formed using the R-algorithm □. This algorithm, however, is not exact: it 
causes systematic errors of order At 2 , where At is the step-size of the Molec- 
ular Dynamics evolution. A careful extrapolation to zero step-size is therefore 
needed to obtain exact results. Nevertheless, it is common practice to forego 
this extrapolation and to perform simulations with a single step-size chosen 
small enough that the expected systematic errors are smaller than the statisti- 
cal ones. We want to point out that there is an alternative to the R-algorithm, 
which gives arbitrarily accurate results without any extrapolation! 

Some time ago, a local algorithm, the so-called "Multiboson algorithm", 
was proposed by LiischerQ, in which the inverse of the fermion matrix is ap- 
proximated by a suitable Chebyshev polynomial. Originally he proposed it for 
two flavors QCD. Borici and de ForcrandH noticed that the determinant of a 
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fcrmion matrix can be written in a manifestly positive way using a polynomial 
approximation, so that one can simulate odd flavors QCD with the multiboson 
method. Indeed, using this method, one flavor QCD was simulated successfully 
I . The same polynomial approximation can be applied for the HMC algorithm 
3, which means that the HMC algorithm also has the possibility of simulating 
odd flavors. Here we give the formulation for simulating odd flavors with the 
HMC algorithm and perform three flavors simulations. Then we compare our 
results with those of the R-algorithm. 

2 Formulation 

2.1 n f =2 

The lattice QCD partition function with nj degenerate quark flavors is given 
by 

Z = J dUdet(D) n f exp(-S gauge ), (1) 

where D is the fermion matrix and in this study we use Wilson fermions. For 
rif = 2, the partition function is 

Z= f dUdet(D) 2 exp(-Sgau g e)- (2) 

In the formulation of the HMC algorithm, det(Z?) 2 is treated as 

dctl? 2 - /d(/) t d0exp(-0 t L> t - 1 L> t 0), (3) 

where the 75 hermiticity of the fermion matrix D, i.e. D = 75-0^75, is used. 

Introducing momenta P conjugate to the link variables U, the partition 
function is rewritten as 

Z= f dUdPcx-p(-H), (4) 

where the Hamiltonian H is defined by 

H = ip 2 + S gauge + tfD^- x D- l <j>. (5) 

This Hamiltonian is used for the Molecular Dynamics (MD) simulation of 
the standard HMC algorithm. Eq.(^) has a computational difficulty in MD 
simulations since one must solve x = D~ l 4> type equations which in general 



take a large amount of computational time for a large fermion matrix and/or 
for a small quark masSj 

Following Luscherel, the inverse of -D can be approximated by a polynomial: 

n 

l/D^P n {D) = Y[(D-Z k ) (6) 

fc=i 

where Zk are roots of the polynomial P n (D): 

2irk 
Z k = 1 - expii — — -). (7) 

n + 1 

The rate of convergence of the approximation depends on the quark mass (See 
Sec.3). 

Replacing D^ 1 in cq.(jg) by P n (D) we obtain an approximate Hamiltonian, 

H n = \p 2 + S gauge + 4tP n {DyP n {D)4>. (8) 

An advantage of using H n is that no solver calculation is required in the MD 
evolution. Instead, one needs n multiplications by the matrix D. Originally 
H n was introduced to reduce computational work. Indeed, iLwas shown that 
H n can provide some gain over the standard HMC algorithm! 

H n does introduce some systematic errors from the polynomial approxi- 
mation. For the n^=even case, however, these errors are easily corrected at 
the Metropolis step by using the exact Hamiltonian of eq. (|jJJ. 

The domain of convergence of P n (D) is bounded by a circle centered at 
(1,0) which goes through the origin. If all eigenvalues of D fall inside this 
domain, P n (D) converges exponentially. Otherwise, P n (D) does not converge, 
which may happen for some exceptional configurations. Our algorithm rejects 
these configurations at the Metropolis step. This domain of convergence can be 
changed by adopting another approximating polynomial. However, the origin 
must be excluded. Together with connectedness and conjugate symmetry of the 
spectrum, this implies that the real negative axis is always excluded from the 
domain of convergence for any polynomial. Configurations with real negative 
Dirac eigenvalues will be rejected by our polynomial algorithm. 

2.2 n f = 1 

In this case, we have to consider det D. det D can not be expressed in a mani- 
festly positive manner using the same treatment of eq. (0) . Thus the standard 
HMC algorithm can not handle n/ = 1 or n/=odd simulations. 



The multiboson algorithm was originally developed for a simulation of 
?i/=2jQCDq. After invention of the multiboson algorithm, Borigi and de For- 
crandtl noticed that a single det D can be treated in a manifestly positive way 
and an nf = 1 multiboson simulation was performed to study thermodynamics 
of n/ = 1 QCDH. 

As before, the inverse of-thc fermion matrix D, using a polynomial of 
degree 2n, is approximated asoD 

2 n 

l/D*l[(D-Z k ), (9) 

fe=i 

where Z k = 1 — exp(i 27rfc/(2n + 1)). Noticing that the Z k s come in complex 
conjugate pairs, eq.(||) is rewritten as 

n 

l/D^H(D-Z k )(D-Z k ). (10) 

Using the 75 hermiticity of the fermion matrix, we find that det(_D — Z k ) = 
det(D — Z k Y . Thus the determinant of D is written as 

det(£>) - det(Tt(£>)T n (D))- 1 , (11) 

where T n (D) = YYk=i(D — Z k ). Using an integral form of the determinant, we 
obtain 

det(D)~ f d(j)U<j)eM-^T 7 \{D)T n (D)(f>). (12) 



The term 0^T r [(D)T n (D)(/) is manifestly positive. Then we may define the 
Hamiltonian of n/ = 1 QCD as 

H = \P 2 + Sgauge + cf>^(D)T n (D)<f>. (13) 

With this Hamiltonian there is no difficulty to perform HMC algorithm. The 
domain of convergence of the approximation eq.(|lj) is the same as for n/ = 2. 
Exceptional configurations for which eigenvalues fall outside this domain will 
be rejected at the Metropolis step. 

2.3 n f = 2 + 1 

The partition function of nf = 2 + 1 QCD is given by 

Z= dU det D 2 det Dexp(-S gauge ), (14) 

1 
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Figure 1: (left): X„ versus degree n. (right):|X n — X exac t\ versus degree n. 



where the notations D and D are introduced to distinguish the two different 
quark masses. Using eq.(0) for detD 2 and eq.rtlj) for detl?, 



det D 2 det D 



expi-^D^D- 1 ^ - <ftT*(D)T n (D)<f>). (15) 



We define rif — 2 + 1 Hamiltonian by 



H = -P 2 + S g + ft&- l D- l 4> + tfTl(D)T n (D)cf>. 



(16) 



Two remarks are in order: (?) the degree n of the approximating poly- 
nomial may be different during the Molecular Dynamics trajectory and for 
the Metropolis acceptance test; the former can be made arbitrarily small and 
tuned for maximum efficiency, while the latter should be taken very large to 
enforce the correct measure; (ii) the two bosonic fields <fi and <j> could be re- 
placed by a single one, with action (j^T^{D)D^~ 1 D~ x T n {D)<f). For simplicity, 
in this exploratory study we use two distinct bosonic fields and a single degree 
for the approximating polynomial. 



3 Convergence 

3.1 n f = 2 

In order to see the rate of convergence of P n (D), we calculate the quan- 
tity X n = $ P^(D)P n (D)4>. In the limit n -> oo, X n goes to X exact = 
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Figure 2: (left): X n versus degree n. (right):|X n — X max \ versus degree n. 

<fi D^ D~ x §. First, we choose X exact — rfi} where r\ is a random gaus- 
sian vector. Then the vector <fi is set to <f> = Dr\. The accuracy of X n is 
measured by the difference between X n and X exact . We use a random gauge 
configuration for this analysis. Figure 1: (left) shows X n versus the degree n 
for different quark masses. Here the same 77 is used for each calculation of X n . 
X n converges to one value as n increases, but at high degree n, X n diverges, 
which can be understood due to the rounding errors of our computer, where 
calculations are performed with 64-bit accuracy. 

Figure l:(right) shows the accuracy of X n by \X n — X exact \. Exponential 
convergence is seen for each quark mass, but the rate of convergence is slow 
for small quark masses. 



3.2 n f = 1 

We do the same analysis as for n/ = 2, but for nf = 1, the value of X exact is not 
known. So we calculate the quantity X n = <ftT^(D)T n (D)4>, where the vector 
is a gaussian random vector, and we use a random gauge configuration. We 
assume that X n goes to a certain value in the limit of n — > 00. 

Figure 2: (left) shows X n as a function of degree n. X n seems to converges 
to a certain value when the degree n increases. At high degree n, X n diverges 



as in the case of n 
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To see the rate of convergence, we calculate \X n — X max \ where X max is 



defined by X„ 



X m , m >> n. Due to the rounding errors, we can not 
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Figure 3: (left): Plaquette of nf = 3 flavor QCD on an 8 2 X 10 X 4 lattice at (3 = 5.0 and 
at re = 0.130 as a function of degree n. (right): Real part of Polyakov loop. 



take very large m. We take a maximum number m where the rounding errors 
still do not appear. Figure 2:(right) shows \X n — X max \ as a function of degree 
n. The dips seen in the figure are just due to the fact that at those points 
X n = X max = X m . The convergence seems to be exponential, but the rate of 
convergence is slow for small quark masses as in the n/ = 2 case. 

4 Simulations 

We perform simulations of three flavors QCD on an 8 2 x 10 x 4 lattice at 
(3 = 5.0 with k = 0.130 and 0.160. We measure the plaquette and Polyakov 
loop varying the degree n and compare them with those from the R-algorithm 
obtained with a step-size At = 0.0113. Figures 3 and 4:(left) show the plaquette 
as a function of n at n = 0.130 and 0.160, respectively. Except for very small 
n, the results from the HMC algorithm agree with those from the R-algorithm 
within statistical errors. Results of the Polyakov loop are shown in figures 3 
and 4: (right). Except for a small discrepancy seen in figure 3, the results from 
the HMC algorithm are in agreement with those from the R-algorithm. Note 
that convergence is not monotonic in n. 



5 Conclusions 

We formulated an odd-flavor HMC algorithm using a polynomial approxima- 
tion. Simulations of three flavors QCD were performed. We found that the 
plaquette values are consistent with those from the R-algorithm at very small 
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Figure 4: (left): Plaquette of nf = 3 flavor QCD on an 8 2 X 10 X 4 lattice at (3 = 5.0 and 
at k = 0.160 as a function of degree n. (right): Real part of Polyakov loop. 



step-size. In principle the HMC algorithm is able to simulate any flavors of 
QCD, with arbitrary accuracy and without extrapolation [as long as all Dirac 
eigenvalues are not real negative]. However the rounding errors should be un- 
der control when we use a large lattice or/and small quark masses where one 
may need a polynomial of high degree n to achieve sufficient approximation. 
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